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Abstract 

A  low  velocity  near  surface  anomaly  delays  the  seismic  wavefield 
and  focuses  wavefronts  passing  through  it.  Such  anomalies  occur  as 
a  result  of  localized  gas  seepages  (“gas  clouds”)  in  the  sedimentary 
column,  for  example.  Use  of  a  migration  velocity  model  not  containing 
the  anomaly  will  produce  images  containing  false  synclinal  structures. 
Application  of  DSO  velocity  inversion  to  a  simple  example  of  this  type 
images  the  velocity  anomaly  and  removes  the  false  structure. 
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Introduction 


Localized  anomalies  in  otherwise  mildly  structured  subsurface  velocity  fields 
pose  an  interesting  challenge  to  velocity  estimation  from  reflection  data. 
Identification  of  anomalies  is  both  simpler  than  the  estimation  of  general 
laterally  heterogeneous  structures  and  more  complex  than  inversion  of  lay¬ 
ered  models.  Anomalies  may  be  fast,  modeling  permafrost  or  salt  intrusions, 
or  slow,  modeling  localized  gas-bearing  rocks.  An  example  of  a  slow  anomaly 
caused  by  the  presence  of  gas  sands  and  its  effect  on  reflection  field  data  ap¬ 
pears  in  Biondi’s  recent  work  [1].  Biondi  uses  a  semblance  measure  based  on 
beam  stacks  to  estimate  the  heterogeneity  and  correctly  image  the  structure 
under  the  sands. 

This  paper  treats  a  simple  synthetic  example  in  which  an  otherwise  lay¬ 
ered  velocity  field  is  altered  by  inclusion  of  a  circular  low- velocity  zone.  The 
physics  of  wave  propagation  used  in  our  experiments  is  2D  constant  density 
primaries-only  acoustics,  which  is  the  simplest  theory  within  which  the  ques¬ 
tions  asked  here  may  be  posed.  The  reflectivity  structure  of  this  model  is  flat. 
However  prestack  migration  (or  linearized  inversion)  of  data  “shot”  over  this 
model,  using  a  layered  velocity  model  not  containing  the  anomaly,  produces 
synclines  and  more  complicated  non-flat  features.  The  object  of  our  excercise 
is:  working  only  with  the  waveform  data,  estimate  the  velocity  sufficiently 
accurately  to  image  the  correct,  flat  underlying  reflector  structure. 

For  this  purpose  we  use  differential  semblance  optimization  (“DSO”),  a 
variant  of  least-squares  inversion.  DSO  has  been  described  at  length  else¬ 
where;  here  we  include  only  a  very  brief  synopsis  of  the  method.  The  DSO 
approach  combines  elements  of  migration  velocity  analysis,  traveltime  to¬ 
mography  and  nonlinear  inversion.  In  particular,  DSO  works  directly  with 
waveform  data,  yet  produces  velocity  estimates  closely  resembling  the  output 
of  traveltime  inversion.  DSO  is  based  on  a  variational  principle.  Optimum 
estimates  are  approximated  by  a  special  iterative  minimization  algorithm. 

For  the  example  constructed  here,  DSO  succeeded  in  removing  almost  all 
of  the  false  structure  in  two  iterations.  The  velocity  anomaly  is  also  imaged. 
The  reconstructed  velocity  shows  the  hallmarks  of  tomographic  inversion: 
smearing  along  the  source/receiver  raypaths,  oscillation  in  the  orthogonal 
direction.  This  is  so  even  though  no  times  are  picked,  no  rays  are  traced  and 
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no  raypath  backprojections  are  performed  (explicitly).  Only  combinations 
of  simulations  and  migrations  are  used  in  the  DSO  process.  Nonetheless 
DSO  produced  a  kinematically  correct  velocity,  of  much  the  same  quality  as 
is  obtained  by  a  functional  traveltime  inversion.  In  this  case  at  least,  DSO 
performed  “tomography  without  picking”. 
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Velocity  inversion  via  DSO 


This  section  gives  a  brief  overview  of  the  2D  primaries  only  constant  density 
acoustic  model  and  differential  semblance  optimization  in  that  context.  For 
a  more  complete  account  see  [6],  [7],  and  [3]. 

The  principal  components  of  the  acoustic  model  discussed  here  are 

•  p(x,z,t)  —  pressure  field 

•  po(x,z,t)  =  reference  pressure  field 

•  v(x,z)  =  velocity  field  (smooth) 

•  r(x,z)  =  reflectivity  field  (~  2 8v/v)  (rough) 

In  addition  we  regard  as  known: 

•  /(0  =  source  time  function  (oscillatory) 

•  i,£  {source  locations} 

•  xT  €  {receiver  locations} 

•  t  €  {time  interval} 

The  source  time  function  f(t )  could  be  included  amongst  the  unknowns 
as  well,  at  the  cost  of  added  computational  expense. 

These  quantities  are  connected  by  the  coupled  wave  equations 

(^)2^2  “  V^)  **(*•’ *•*)  =  fi*)S(X  ~  X°) 

fid2  \ 

UAxydt*  ~  V*J  =  r(x)V2po(xs,x,t) 
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plus  appropriate  initial  and  boundary  conditions.  The  first  equation  is 
the  usual  acoustic  wave  equation  for  radiation  from  an  isotropic  point  source. 
The  second  is  the  perturbational  equation  for  pressure  fluctuations  p(x ,  z,t ) 
due  to  the  velocity  perturbation  r  =  8v/v,  i.e.  the  primary  reflection  or 
singly  scattered  field. 

Reflection  data  is  predicted  by  the  forward  map : 

F[u,r]  =  {p(xs,:rr,t)} 

For  simplicity  receiver  arrays  are  simplified  to  points  {xr}  as  well. 

The  reflection  data  inversion  problem  may  be  stated: 

given.  Fdata  —  {Pdata(^'S)  2-rj 
find:  velocity  v;  reflectivity  r 
so  that 


F[v,  r]  ~  Fdata 


As  have  many  others,  we  formulate  this  task  as  a  best-fit  problem,  via  the 
least  squares  principle  (see  e.g.  [9]): 

choose  v,  r  to  minimize 
||F[v,r]  -  Fdata II2 

(||...||  =  L 2  norm  =  mean  square  fit  error) 

The  nature  of  this  variational  problem  follows  from  key  properties  of  the 
forward  map  ,F[t;,r]  =  {p(a:s,  xr,  f)},  which  is 

•  linear  in  r 

•  very  nonlinear  in  v 

Consequently  the  mean  square  fit  error  ||F[u,  r]  —  Fdata\\2  is 

•  quadratic  in  r,  but 
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•  highly  nonquadratic  (nonconvex)  in  v 

Critical  consequences  for  the  least  squares  principle  are: 

•  estimation  of  r  (oscillatory)  given  v  (smooth)  is  relatively  easy:  there 
are  many  algorithms,  and  all  work  relatively  well; 

•  estimation  of  v  and  r  simultaneously  is  very  difficult; 

•  least-squares  inversion  for  v  appears  to  requires  nonconvex  (global) 
optimization  (Monte-Carlo,  simulated  annealing...). 

Many  other  velocity  indicators,  for  example  the  widely  used  NMO-based 
velocity  spectrum  [8]  and  Biondi’s  beam  stack  semblance  [1],  share  the  highly 
nonquadratic  nature  of  the  least  squares  principle  for  velocities.  These  are 
not  the  only  variational  principles  one  can  imagine  for  determination  of  v 
and  r  however.  For  reasons  discussed  extensively  in  the  references,  we  have 
pursued  a  different  variational  principle,  differential  semblance  optimization , 
or  DSO.  Its  key  properties  are: 

•  it  is  a  modification  of  the  least-squares  principle,  and  can  in  fact  be 
regarded  as  an  infeasible  point  method  for  finding  the  least  squares 
solution; 

•  it  regards  the  model  for  each  experiment  (x»)  as  independent: 

r  =  r(xs,x)  (but  v  —  u(:r)) 

•  since  in  fact  the  model  should  not  be  xs- dependent  (infeasible  -  “there 
is  only  one  earth!”),  a  penalty  for  infeasibility  is  added  to  the  mean 
square  fit  error  to  produce  the  augmented  objective  function 

^{|F[v,r]  -  Fdata |2  +  a2\drldxs\2} 

•  Minimization  of  this  objective  is  accomplished  in  two  stages:  minimize 
first  over  r  (in  which  the  objective  is  quadratic!)  to  produce  a  function 
of  v : 

J<j[v]  =  minimumr(X3tX)^{\F[v,r\  -  Fdata\ 2  -f  cr2\dr / dxs\2} 
then  minimize  Jc  over  v. 
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•  Note  that  computing  Ja  requires  solution  of  “inner”  inverse  problem 
for  r  given  v. 

The  first  term  of  the  differential  semblance  objective  forces  the  model  to 
fit  each  shot  gather  independently.  The  second  term  measures  the  semblance 
of  neighboring  reflectivity  inversions  in  a  differential  sense;  in  practice  a  di¬ 
vided  difference  is  used  in  place  of  the  derivative  with  respect  to  xs.  Since 
neighboring  reflectivity  gathers  are  compared,  reflectivity  gathers  are  similar 
even  when  v  is  substantially  wrong.  Thus  the  differential  semblance  measure 
changes  much  more  smoothly  than  do  global  measures  such  as  stack  power. 
Inverted  reflectivities  are  compared,  rather  than  migrated  images  with  un¬ 
controlled  amplitudes.  Therefore  simple  differences  are  sufficiently  robust  to 
measure  semblance.  This  could  not  be  the  case  if  migrated  images  were  used 
instead,  as  in  migration  velocity  analysis. 

Analysis  of  and  numerical  experiments  with  DSO  show  that 

•  Ja  is  smooth  and  convex  over  a  large  domain  in  model  space  for  small 
a2 

•  if  data  noise  is  small,  the  global  minimizer  of  DSO  is  close  to  the  global 
solution  of  least-squares  problem 

•  The  DSO  Hessian  is  closely  related  to  the  Hessian  of  the  travel  time 
tomography  objective  function  (this  observation  is  due  to  Hua  Song). 

Our  implementation  includes  special  devices  to  ensure  accurate  calcula¬ 
tion  of  the  gradient  of  J0.  This  calculation  is  the  basis  of  a  nonlinear  con¬ 
jugate  gradient  algorithm  using  the  Polak-Ribiere  step  update  formula  ([2]) 
for  minimizing  the  DSO  objective  function.  More  discussion  and  a  detailed 
description  of  the  implementation  of  DSO  used  here  is  found  in  [3]. 
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The  Model 


The  velocity  model  used  in  our  experiments  was  based  on  a  layered  back¬ 
ground  vstrat ,  profiled  in  Figure  1,  with  velocities  ranging  from  1.5  m/ms 
near  the  surface  to  1.9  m/ms  at  1200  m  depth.  This  background  model  was 
depressed  by  roughly  12%  in  a  circular  region  of  diameter  400  m  centered  500 
m  below  the  surface,  to  produce  the  target  velocity  model  vexact  field  shown 
in  Figure  2.  The  reflectivity  field  rexact  used  with  vexact  to  generate  the  data 
is  horizontally  stratified,  with  profile  shown  in  Figure  3. 

A  small  reflection  survey  of  30  shots  over  this  model  was  simulated  using 
a  finite  difference  method  of  fourth  order  accuracy  in  space  and  second  order 
accuracy  in  time.  The  shot  spacing  was  100  m  and  shot  depth  was  8  m.  The 
first  shot  was  located  roughly  1400  m  to  the  west  (left)  of  the  anomaly  center, 
the  last  an  equal  distance  to  the  east  (right).  The  source  was  punctual,  as 
described  in  the  last  section.  Thirty  four  point  receivers  spaced  50  m  apart 
at  12  m  were  distributed  between  150  m  and  1800  m  offset  to  the  east  (right) 
of  the  shot  points.  Sample  intervals  in  space  for  the  simulation  were  16  m  in 
both  directions.  The  sample  interval  in  time  was  4  ms.  The  record  length 
was  2000  ms.  The  source,  regarded  as  known  in  these  experiments,  was  a 
Ricker  wavelet  of  peak  frequency  15  Hz. 

Selected  shot  gathers  are  displayed  in  Figure  4.  The  moveout  distortion 
due  to  the  velocity  anomaly  is  clearly  visible  and  moves  across  the  gather 
as  the  shot  position  varies  from  west  to  east  over  the  anomaly.  A  stack 
of  the  30  shot-dependent  reflectivities  obtained  by  inversion  at  the  layered 
background  velocity  vatrat  of  Figure  1  is  displayed  as  Figure  5.  This  stack  may 
be  viewed  as  a  prestack  migration  image,  and  is  very  similar  kinematically  to 
the  prestack  least  squares  inversion  diplayed  in  Figure  6.  Both  images  show 
the  expected  pull-down  below  the  anomaly,  which  has  become  a  sinusoidal 
undulation  at  the  strong  reflector  at  1300  m  due  to  the  formation  of  a  caustic. 
A  similar  effect  in  the  migration  of  a  field  data  set  may  be  observed  in  Figures 
4,  6,  and  7  of  [1]. 
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DSO  inversion 


Several  parameter  choices  are  required  by  DSO  at  this  stage  in  its  develop¬ 
ment.  It  seems  likely  that  many  of  these  choices  could  be  made  automatically, 
but  at  present  all  are  made  by  trial  and  error.  Two  of  the  most  important 
of  these  paramters  are  the  differential  semblance  weight  a  and  the  length 
scale  of  the  allowed  velocity  models.  After  some  experimentation  we  chose 
o'  =  0.01.  The  length  scale  (i.e.  degree  of  smoothness)  of  the  velocity  model 
is  controlled  by  choice  of  norm  in  velocity  space,  rather  than  by  parsimonious 
parameterization  (e.g.  [1])  to  avoid  bias.  The  norm/smoothing  mechanism 
is  described  in  [3].  We  chose  norm  parameters  to  suppress  wavelengths  in 
the  velocity  shorter  than  roughly  400  m  in  these  experiments. 

Figure  7  displays  some  common  depth  point  gathers  (or  common  reflec¬ 
tion  point  gathers  or  coherency  panels)  for  the  reflectivity  estimate  at  vstTat- 
Considerable  moveout  is  visible;  just  as  in  migration  velocity  analysis,  this 
indicates  a  need  to  update  the  velocity  model  (see  e.g.  [10]  for  many  examples 
of  this  idea). 

The  gradient  of  the  reduced  DSO  objective  function  at  the  layered  back¬ 
ground  velocity: 


gradJ^Vstrat) 

is  shown  in  Figure  8.  The  vertical  smearing  and  side  lobes  are  characteristic 
of  limited  aperture  tomography;  compare  for  example  [5],  upper  left  part  of 
Figure  6  which  has  similar  geometry.  Nontheless  the  location  and  extent  of 
the  anomaly  is  already  well  defined.  Note  that  no  a  priori  information  about 
either  the  location  or  size  of  the  anomaly  has  been  provided. 

Two  steps  of  the  nonlinear  conjugate  gradient  algorithm  produced  the 
velocity  model  r>2  displayed  in  Figure  9.  The  low  velocity  zone  is  correctly 
identified,  modulo  the  vertical  smearing  noted  earlier.  The  reflectivity  CDP 
gathers  for  this  model  (Figure  10)  show  that  most  of  the  residual  moveout 
has  been  removed.  Moreover  the  stacked  reflectivities  (Figure  11)  are  much 
flatter;  the  maximum  deflection  in  the  center  of  the  synclinal  structures  has 
decreased  from  about  50  m  in  Figure  6  or  7  to  less  than  10  m  in  Figure  9, 
which  is  close  to  optimal  given  the  bandwidth  of  the  data.  It  seems  likely 
that  further  iterations  could  remove  at  least  some  of  the  very  mild  remaining 
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structure,  which  is  now  roughly  10%  of  a  wavelength  in  vertical  extent. 


Conclusion 

This  synthetic  example,  while  extremely  simple,  shows  that  DSO  is  capable 
of  resolving  lateral  velocity  heterogeneities.  As  we  have  shown  elsehwere  ([3]) 
that  DSO  effectively  estimates  layered  structures  beginning  with  quite  inac¬ 
curate  starting  models,  we  might  hope  that  DSO  will  provide  kinematically 
correct  estimates  of  rather  arbitrarily  heterogeneous  velocity  fields. 

A  major  drawback  of  the  current  implementation  is  its  computational 
expense:  the  two  iterations  used  here  required  roughly  140  simulations  and 
migrations  of  the  entire  data  set,  and  cost  about  10  Cray  Y-MP  CPU  hours. 
This  cost,  and  even  the  projected  cost  of  more  accurate  inversions  with  higher 
frequency  content  and  larger  models,  is  far  smaller  than  the  projected  cost  of 
stochastic  least  squares  inversion  as  proposed  for  instance  by  [4].  Nonetheless 
it  is  larger  than  one  would  like,  and  will  hamper  experimentation  for  some 
time  to  come  unless  reduced. 

In  the  experiments  reported  here,  more  than  95%  of  the  CPU  cycles  were 
spent  in  finite  difference  simulation  and  migration,  which  runs  at  about  250 
Mflops  on  the  Y-MP.  This  is  roughly  70%  of  the  peak  floating  point  through¬ 
put  of  the  Y-MP  CPU;  somewhat  higher  Mflops  ratings  are  achieved  with 
larger  models,  though  the  total  CPU  time  goes  up  considerably.  Therefore 
the  speed  of  the  finite  difference  portion  of  the  code  is  essentially  hardware- 
limited,  and  we  must  look  elsewhere  for  significant  software- driven  speedup. 
We  are  experimenting  with  improvements  in  the  optimization  algorithm, 
preconditioning  of  the  inner  (reflectivity)  inversion,  and  parallelization  over 
shots  on  both  distributed  and  shared  memory  architectures,  amongst  other 
devices  all  aimed  at  a  considerably  faster  implementation. 
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Figure  Captions 


1.  Profile  of  stratified  background  velocity  model  vstrat- 

2.  Velocity  model  vexact  used  to  generate  data.  Near-surface  velocity  is  1.5 
m/ms;  lowest  velocity  in  upper  part  of  anomaly  is  roughly  1.4  m/ms. 
High  velocity  below  1100  m  is  1.9  m/ms.  Horizontal  tick  marks  are 
trace  number;  traces  are  spaced  16  m  apart. 

3.  Profile  of  stratified  reflectivity  model  rexact- 

4.  Selected  shot  gathers.  Vertical  scale  is  ms;  trace  offset  interval  is  50 
m,  far  offset  is  1800  m.  A  mute  has  been  applied  to  remove  direct  and 
refracted  energy. 

5.  Stack  of  inverted  reflectivities  at  the  stratified  velocity  model  of  Figure 
1.  Scales  are  same  as  Figure  2.  Note  distinct  pulldown,  becoming 
sinusoidal  below  a  coustic  at  the  bottom  of  the  display. 

6.  Result  of  least  squares  inversion  (also  performed  with  the  DSO  software 
-  it  does  both!).  Compare  Figure  5. 

7.  CDP  gathers  of  reflectivity  inverted  at  stratified  model,  from  1000  m 
west  to  1000  m  east  of  the  anomaly.  CDP  separation  is  200  m.  Note 
the  considerable  moveout  over  the  center  section,  which  is  affected  by 
the  anomaly.  This  is  a  clear  indication  of  velocity  error. 

8.  Gradient  of  the  DSO  objective  function  Ja  at  the  stratifiec  background 
velocity  v3trat  (Figure  1).  Comparing  with  Figure  2,  see  that  anomaly  is 
already  located  rather  precisely.  Typical  tomographic  raypath  smear¬ 
ing  and  side  lobes  are  present.  Scales  same  as  Figure  2. 

9.  Velocity  estimate  after  2  nonlinear  conjugate  gradient  steps.  Compare 
Figure  2. 

10.  CDP  gathers  of  reflectivity  inverted  at  model  of  Figure  9.  Compare 
Figure  7;  note  that  gathers  are  essentially  flat. 

11.  Stacked  reflectivities  inverted  at  model  of  Figure  9.  Compare  Figures 
5  and  6;  note  that  false  structure  is  almost  entirely  removed. 
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